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Abstract — Unsteady flowfields of a two-dimensional oscillating airfoil are calculated using an implicit, 
finite-difference, Navier- Stokes numerical scheme, Five widely used turbulence models are used with the 
numerical scheme to assess the accuracy and suitability of the models for simulating the retreating blade 
stall of helicopter rotor in forward flight. Three unsteady flow conditions corresponding to an essentially 
attached flow, light-stall, and deep-stall cases of an oscillating NACA 0015 wing experiment were chosen 
as test cases for computations. Results of unsteady airloads hysteresis curves, harmonics of unsteady 
pressures, and instantaneous flowfield patterns are presented. Some effects of grid density, time-step size, 
and numerical dissipation on the unsteady solutions relevant to the evaluation of turbulence models are 
examined. Comparison of unsteady airloads with experimental data show that all models tested are 
deficient in some sense and no single model predicts airloads consistently and in agreement with 
experiment for the three flow regimes. The chief findings are that the simple algebraic model based on 
the renormalization group theory (RNG) offers some improvement over the Baldwin Lomax model in 
all flow regimes with nearly same computational cost. The one-equation models provide significant 
improvement over the algebraic and the half-equation models but have their own limitations. The 
Baldwin-Barth model overpredicts separation and underpredicts reattachment. In contrast, the 
Spalart- Allmaras model underpredicts separation and overpredicts reattachment. 


INTRODUCTION 

The term ‘'dynamic stair’ is usually referred to the unsteady separation and stall phenomena of 
aerodynamic bodies or lifting surfaces that are forced to execute time-dependent (unsteady) motion, 
oscillatory or otherwise. It is a complex fluid dynamic phenomenon of practical importance and 
occurs on maneuvering flight vehicles, retreating blades of helicopter rotors, wind turbine blades, 
and compressor cascades. The dynamic stall phenomena often lead to the initiation of stall flutter. 
As summarized in extensive reviews by McCroskey [1, 2] and Carr [3], the majority of the work 
on this fundamental fluid dynamic problem is devoted to the case of airfoils oscillating with 
moderate amplitude in a uniform freestream. When the airfoil reaches fairly high angles of attack 
during the oscillatory cycle, past the static stall angle limit, the generated unsteady flowfield is 
characterized by massive separation and large-scale vortical structures. One important difference 
between this flowfield structure and that generated by the static stall is the large hysteresis in the 
unsteady separation and reattachment process. The maximum values of lift, drag, and pitching- 
moment coefficients can greatly exceed their static counterparts, and not even the qualitative 
behavior of these can be reproduced by neglecting the unsteady motion of the body surface. 

One of the reasons why the flowfield associated with dynamic stall is more difficult to analyze 
than the static stall is of its dependence on a much larger number of parameters. The important 
ones are the airfoil shape, Mach number, reduced frequency or pitch rate, amplitude of oscillations, 
type of motion (ramp or oscillatory), Reynolds number, three-dimensional effects, and wind tunnel 
effects. To date most of the research in this area has been performed for the simpler model problems 
of two-dimensional oscillating airfoils. Most of what is understood about the characteristics and 
various regimes of dynamic stall has essentially come from experimental observations, which are 
mostly two-dimensional. Attempts to calculate the quantitative effects of dynamic stall abound in 
the literature [3-1 1]. The purely laminar case is assumed to have been solved [5, 6], although recent 
studies [12] show small-scale details of possible importance. However, the laminar calculations have 
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not been validated with experiments (for the lack of availability of purely laminar data). The flows 
with turbulent boundary layer have not yet been successfully solved. 

The weak link in the computational and theoretical methods for an accurate simulation of these 
unsteady flowfields is the turbulence modeling. Of course, the transitional nature of the boundary 
layer is almost always neglected; instead the flow is approximated to be either completely laminar 
or completely turbulent on the airfoil surface. Such an approximation may not be correct if the 
flowfield is dominated by leading-edge separation. In any case, a reasonably good turbulence model 
must be used to accurately calculate the nonequilibrium nature of the separated turbulent boundary 
layer and the associated unsteady time-lag features. Simple eddy viscosity models, such as the 
Baldwin-Lomax model [13], have been found to be inadequate. The objective of this investigation 
is to identify a reasonably accurate and robust turbulence model. Several turbulence models that 
are in use in most Computational Fluid Dynamics (CFD) codes for the computation of steady-state 
flows are considered for evaluation. It should be noted that Refs [7-10] have considered a similar 
exercise. 

The five turbulence models considered are the Baldwin-Lomax (B-L) algebraic model [13], the 
Renormalization Group theory (RNG) based algebraic model [14], the half-equation Johnson-King 
(J-K) model [15, 16], and the one-equation models of Baldwin-Barth (B-B) [17] and 
Spalart-Allmaras (S— A) [18]. The performance of these models is evaluated for accuracy and 
robustness by using them to calculate the unsteady, two-dimensional, viscous, flowfields of an 
oscillating NACA 0015 airfoil. The parameters of grid size, numerical dissipation, and time-step 
size was varied to arrive at a set of these parameters for accuracy and robustness. The accuracy 
is validated by comparison with oscillating wing experimental data [19] measured at the U.S. Army 
Aerofiightdynamics Directorate at the NASA Ames Research Center. The eventual objective of this 
study is to identify a turbulence model that calculates accurately the flow physics of unsteady 
boundary layer, its separation and reattachment process, and that is appropriate for modeling the 
three-dimensional retreating blade stall of a helicopter rotor in forward flight. 


GOVERNING EQUATIONS 

The governing equations considered are the Reynolds-averaged, two-dimensional, 
Navier-Stokes equations in strong conservation-law form. These can be written in a generalized 
body-conforming curvilinear coordinate system (£,>/, t) as follows [20]: 

d,Q + d ( £ +d,F = 1.(8^ + s,S) ( 1 ) 

where £, = £(jc, y, t), rj = t](x,y, t ), and x = t. Here (x, y, t) is the inertial coordinate system; x is 
in the streamwise direction and y is normal to it. Also, Q is the vector of conserved flow variables; 
E and F are the in viscid flux vectors. These are given by 
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The vectors /? and S are the viscous stress vectors in the £ and q directions, respectively. The viscous 
terms may be retained in both directions to resolve massive separation and these are considered 
in the thin layer approximation. For example, the vector in the -direction is written as 
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Here 


m x = t)l + ti 2 y 


m 2 = rixU n + riyV n 


m,= 



(i u 2 + v 2 )/2 + 


1 

Pr(y - 1) 



w 4 = rj x u + t} y v 

Also, U and V are the contravariant velocity components. The terms £ y , rj y , and r\, are 
the coordinate transformation metrics and J is the Jacobian of the transformation [20], 

In the above equations, all geometrical dimensions are normalized with the airfod chord c, the 
Cartesian velocity components u and v are scaled by the freestream sound speed a x , and the time 
t is normalized as tcja x \ p is the static pressure normalized by p x a 2 x ; p is the density normalized 
by free-stream density p x ', e is the total energy per unit volume normalized with p x a x ., a is the 
speed of sound; Re is the Reynolds number; Pr is the Prandtl number; y is the ratio of specific heats; 
and p is the viscosity coefficient normalized by its free-stream value. The pressure is related to the 
density and total energy through the equation of state for an ideal gas given by, 

p = (y - 1 ) [e - p(u 2 + v 2 )/2] (2) 

Finally, the rotational speed of the airfoil, «, is obtained from the type of motion prescribed 
as co = da/d t. The reduced frequency parameter is defined as either k = dc/2U x , or nfcjU x where 
/ is the frequency of oscillation and U x is the free-stream velocity. Then 
w = a, (2kU x /c)cos[(2kU x /c)t] for «(?) = «« + *i sin(cof) where <x 0 is the mean angle of oscillation 
and a, is the amplitude of pitch. 


NUMERICAL SCHEME 


The governing equations are solved using a finite-difference, implicit, approximately factored 
Beam-Warming numerical algorithm [20]. The viscous terms are treated implicitly and they are 
retained in both and q -directions and in a thin layer approximation form. The approximately 
factored algorithm is given by 


{/ + h[8 f A", - Rc 'd.M/, + (Amp ),]}''! 7 + - R e~ l 5 n % + (Amp)J}'’ A 0C, 

= - a/ { 0? y - 0 h + [» e K + - VtK + (3) 


Here d is the central difference operator and h is the time-step factor that determines whether the 
algorithm is first- or second-order time accurate. The time mdex is denoted by n and 
&Qn = &,t 1 — &j)- The explicit inviscid fluxes are given by E u , F u and R^, 3,j are the explicit 
viscous fluxes. The quantities. A, S, M, and N in equation (3) are flux Jacobian matrices obtained 
from the linearization of the left hand side and A mp an d are the implicit and explicit dissipation 
terms. A Jameson-type [21], blended second and fourth order numerical dissipation, based on the 
computed pressure field, is used to suppress high frequency numerical oscillations. For subsonic 
shock free solutions, only the fourth-order dissipation is used, while for transonic solutions the 
second-order dissipation is activated in the vicinity of shocks where the pressure jump has steep 
gradients. Both implicit and explicit dissipation terms are scaled by the spectral radius. For the 
accuracy of calculated solutions, the added dissipation coefficients are kept as small as possible. 

The errors introduced by the linearization and approximate factorization of the left hand side 
of the numerical algorithm may be minimized by performing Newton subiterations at each 
time-step during the unsteady calculations. The approximation to Q" +l at each subderation is the 
quantity Q p . When p >2, during a given level of subiteration to convergence, Q p =Q n+ \ but when 
p = 1 and no subiterations are performed Q” = Q", and Q p+ ' = Q n+ '. In the present study, the 
numerical experiments have demonstrated that because of the small time-steps used, the Newton 
subiterations are not required. It was also found that the two time-level numerical scheme does 

not increase the accuracy of the unsteady calculations. 

In the normal practice of the thin layer approximation for viscous terms, only the terms in the 
normal direction (normal to the solid surface) are retained because of the large flow gradients in 
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that direction. Retaining the viscous terms in both directions and in their thin layer form was found 
to be slightly beneficial for the deep-stall cases, perhaps because of massive separation. For the 
light-stall case, however, the calculations performed retaining first the viscous terms for the 
^-direction alone and then the terms for both £- and >7 -directions showed very little difference 
between the solutions. Therefore, the light-stall calculations are performed by retaining the viscous 
terms only in 17 -direction. The computational cost was not significantly increased for doing this 
for both the directions. Also, the numerical experiments have demonstrated that implicit treatment 
of the streamwise viscous term, M of equation (3), does not contribute to the accuracy of the 
solution but results in increased computational cost. Therefore, implicit treatment of the streamwise 
viscous term is not used in the present study. 

Body-fitted C-type computational grids are used in all calculations. These are generated using 
a hyperbolic grid generator. The grids are clustered at the body surface in the normal direction, 
leading edge, and trailing edge regions. The spacing of the first grid point at the surface in the 
normal direction is 0.00002 chord and this translates to a = 0(1) at the mid-chord of the airfoil. 
The grid boundaries are located at 1 5 chords in all directions. Although most of the results 
presented are generated using one grid of size 361 x 71 with 271 points on the airfoil in the 
chord wise direction and the remaining in the wake region, three other grids of size 181 x 71, 

671 x 71, and 361 x 141 have also been used to study the effect of grid size on the solution accuracy! 

In all cases, the grids are designed to have at least 25 points in the boundary layer. 

Boundary conditions are updated explicitly. For subsonic inflow-outflow, the flow variables at 
the boundaries are evaluated using one-dimensional Riemann invariant extrapolation. On the 
airfoil surface the velocities are set equal to the surface velocity. The density and pressure are 
evaluated from extrapolation from its interior neighboring grid point. For C-type grids used in this 
study averaging of the flow variables at the wake cut is used. 

TURBULENCE MODELS 

All flows computed in this study are assumed to have fully turbulent boundary layer on both 
upper and lower surfaces of the airfoil by neglecting the laminar and transitional boundary layer. 
In the experiments of Piziali [19], with which all present calculations are compared, the boundary 
layer is tripped in the leading-edge region. (The experimentalist has performed several checks to 
assertain that the boundary layer downstream of the trip is fully turbulent.) For turbulent viscous 
flows, the non-dimensional viscosity /i in the viscous flux vectors is calculated as the sum total of 
the laminar and turbulent viscosity. It is the determination of this turbulent viscosity that is of 
special significance and the focal point of this study. As mentioned before, five different turbulence 
models are used for calculating turbulent eddy viscosity and the unsteady flowfield of an oscillating 
NACA 0015 airfoil. The results are used to evaluate their performance. The details of the 
formulation of these models are presented elsewhere [13—18], However, the following paragraphs 
describe briefly the salient features of these models and the specific versions used in the present 
investigation. 

Baldwin-Lomax (B-L) model 

This is a two-layer, zero-equation (algebraic) model. It is patterned after Cebeci-Smith model 
[22] and introduces a modification that eliminates the need to search for the edge of the boundary 
layer to determine the length scale. It is the most commonly used turbulence model available in 
most of the CFD codes. Its strength and weakness are well known in CFD community; it predicts 
accurately the steady flows with little or no separation and performs poorly if there is large 
separation, either shock-induced or otherwise. It uses an inner and outer layer formulation for 
determining the turbulent viscosity with a smooth transition that spreads over the two regions. It 
uses a classical mixing-length hypothesis for the inner layer with a van Driest damping function 
to force the eddy viscosity at the wall to zero. In the outer layer, the length scale is fixed by the 
location where the product of distance from the solid wall and vorticity reaches a maximum in the 
boundary layer. The Klebanoffs intermittency factor is used to drive the eddy viscosity to zero 
in the outer flow away from the wall. Some of the constants of the theory are determined by 
correlating with experimental data. The details of the theory are described in Ref. [13], 
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RNG model 

Another algebraic eddy viscosity model was proposed recently, for the closure of the Reynolds 
averaged Navier-Stokes equations, based on the Renormalization Group (RNG) theory of 
turbulence [14]. The algebraic model, although free from uncertainties related to the experimental 
determination of empirical modeling constants, still requires specification of an integral length-scale 
of turbulence, similar to the B-L model, which reduces the generality of the model. In this model 
the integral scale is assumed to be proportional to the boundary layer thickness 3, and the eddy 
viscosity is obtained as in Ref. [14] from 


V 



(4) 


where v = v t + v„ the subscripts t and 1 refer to the turbulent and laminar components, respectively. 
Also, 6 is the boundary layer thickness and it is determined from <5 = 1.2 y l/2 , where y il2 * s th e 
normal distance from the wall at which the vorticity function F(y) (see Ref. [13]) attains its half 
amplitude [16]. H is the Heaviside step function and <f> is the dissipation function. The parameter 
a = 0.0192 corresponding to the von Karman constant k = 0.372, and the parameter C c = 75. The 
turbulent eddy viscosity is then obtained by solving equation (4) at every point in the computational 
domain. In estimating the eddy viscosity with this model in this study, the model is applied only 
to the suction side of the airfoil (upper surface) while the pressure side (having attached flow) and 
wake regions are computed with the B— L model. Application of the model to both upper and lower 
surfaces essentially gave the same results as the one obtained by applying for only upper surface; 
so the latter was used. 


Johnson-King (J-K) model 

The above two models, viz., the B-L and RNG models, are termed equilibrium models meaning 
that the eddy viscosity instantaneously adjusts to the local flow without any history effects. The 
next three models presented are called non-equilibrium models in which the calculated eddy 
viscosity accounts for the upstream history of the flow. 

From the time Johnson and King first introduced their half-equation turbulence model [15], there 
have been several modifications to improve their original model for separated flows [9], In the 
present application to unsteady flows, the original version of this model is used [15] and is briefly 
described in the following paragraphs. 

The Johnson King model [15] takes into account the convection and diffusion effects on the 
Reynolds shear stress —u'w' in the streamwise direction. The eddy viscosity is given by 

*-\[i -«■{-$] <5 > 

where v t , describe the eddy viscosity variation in the inner and outer part of the boundary layer. 
The inner eddy viscosity is computed as 

v,. = D 2 KyJ{-u'w’) m i X 

£> = 1 — e~ WA+) ( 6 ) 

where the constant A + = 15. The outer eddy viscosity is given by 

= <r(jt)[0.0168£/ c <5*y] (7) 

where 3* is the boundary layer displacement thickness, y is the KlebanofTs intermittency function 
given by y = [1 + 5.5(>\/c)) 6 ]' , » and a(x) is obt ained fr om the solution of an ordinary differential 
equation which describes the development of — u'w'\ miX along the path of the maximum shear 
stress. The effects of convection and diffusion on the Reynolds shear stress development are 
accounted from the solution of the following ordinary differential equation 
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Here C dir and a, are modeling constants ii m is the maximum average mean velocity and 

g=[-u'w' n Y m , and g eq = 2 

Also, L m is the dissipation length evaluated as 

L m = 0.40y for y m /5 <0.225 

L m = 0.09(5 for yJS >0.225 (9) 

The boundary layer thickness, <5, is determined in the same way as explained in the discussion of 
the RNG model. The equilibrium shear stress g eq in equation (8) is determined from the following 
equilibrium eddy viscosity distribution 



Vi ~ D ~ Ky \/(-iv'w')m, cq 

v lo ,eq = 0.0168C/ e <5*y (10) 

where U e is the velocity at the edge of the boundary layer. 

An implicit Euler method is used for the numerical solution of equation (8), and the maximum 
shear stress at each iteration level is updated as follows 

a(xf*' = (U) 

l \j 

It should be noted that the unsteady term is neglected in the above formulation. Solutions with 
the Johnson-King turbulence model are obtained as follows. First a convergent solution using the 
Baldwin Lomax turbulence model for the entire flowfield is obtained. Then the Johnson-King 
model is applied only to the upper surface of the airfoil as using it for both the surfaces did not 
change the results. To initiate the solution a(x) in equation (7) is set unity and it is allowed to 
change according to equation (11). It should be noted that the Johnson-King model reduces to 
the Cebeci-Smith model [22] when cr(.v) is identically equal to one. 


One - equa t ion models 

The B-L and the RNG models are equilibrium models, in which the production is identically 
equal to the dissipation. The J K model is an improvement over the equilibrium turbulence models 
because it accounts for the evolution of the maximum shear stress through the solution of an 
ordinary differential equation (ODE). It, therefore, attempts to calculate the non-equilibrium 
turbulent boundary layer. I he validity of the models discussed thus far is limited and questionable 
when applied to a flow environment consisting of unsteady separated flow with multiple shear 
layers. 

Recently, several one-equation models have been developed for use in place of these lower-order 
turbulence models. In the present investigation two such models are considered for investigation. 
These are the Baldwin Barth [17] and Spalart-Allmaras [18] models. The primary advantage of 
these models is that they do not require the evaluation of flow-dependent length scales, such as 
the boundary layer thickness. The validity of these models for steady flows has been demonstrated, 
but only in a limited sense. In the present investigation these models are tested for several unsteady 
attached and separated flows over oscillating airfoil. 


Baldwin-Barth (B B) model 

This one-equation model [17] is derived from the simplified form of the k -c model equations. 
It solves a partial differential equation for the modified turbulent Reynolds number vR T from 

D(vRj) t = / Y \ ^ I 

Dl = ta/j - <<, k/v/^ + l V + -± lV 2 (vR T ) - -(Vv,) • V(vR t ) (12) 

This equation solves for the field quantity R T = k 2 /vc = R T / 3 (R T ), named turbulent Reynolds 
number. The turbulent Reynolds number is related to the eddy viscosity as 



Evaluation of turbulence models 


839 


V, = vcj„ R t = vcjji R t ( 1 3) 

where 

- = (c i2 ~C t ,) s /cjK 2 

<*c 

v, = c lx (vRj)D ] D 2 
= P v . 


U = D,D, 

D : = 1 -exp(->' + M + ) 
D , = 1 - exp (-y + /At) 



Here y + = w t >’/v and w t is the skin friction velocity. The constants used for the B-B model are the 
same as in their original paper [17] and are given by: 

k = 0.41, c £| = 1.2, c, 2 = 2.0 £■„ = 0.09, A + = 26, /< 2 + = 10. 

This model is applied to the entire fiowfield to compute the eddy viscosity. 


Spalart A llmar as (S-A ) model 

The second one-equation model used in the present investigation is the Spalart-Allmaras (S-A) 
model [18]. This model requires the solution of a transport partial differential equation for the 
turbulent eddy viscosity. This equation was constructed using empirical criteria and arguments 
from dimensional analysis. It has many similarities with the B-B model. The S-A model also has 
a provision for transition onset at any specified location, although the present investigation treats 
the boundary layer as turbulent on the entire surface. The eddy viscosity is obtained from the 
solution of the following partial differential equation. 

^ = c bl ( 1 — f i2 )Sv + V '((v+ v)Fv) + c b2 (Pv) 2 ] - +/,i^ 2 (14) 

where 

fa = c x i g, exp^ - e l2 [d 2 + g] </?]^ 

/ l2 = c L , exp(-c t4 T 2 ) 

g t = min(0.1, AU/ai^Ax) (13) 

where X = v/v and co, r is used here to denote the vorticity on the wall at the boundary layer trip 
point. The constants of this model have been chosen the same as in the original reference [18], and 
the transition location was set at the airfoil leading edge. The model constants are. 


c bl = 0.1355 c b2 = 0.622 o = 2/3 = c bl /K 2 + (1 + c b2 )/<r c -2 = 2.0 


c„. 3 = 0.3 k = 0.41 f,=g 


1+c 6 b3 v ' 6 


f 6 + c»y 


g = r + c„ 2 (r 6 - r), r = ^~Tj i c n = 1-0 c, 2 -2.0 c, 3 -l.l c l4 - 2.0 
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The computational costs vary widely for these models. The relative costs for these five turbulence 
models is discussed at the end of Results Section. 

RESULTS AND DISCUSSION 

This paper describes an attempt to evaluate five widely used turbulence models in calculating 
the unsteady, two-dimensional flowfields of an oscillating NACA 0015 wing spanning a wind tunnel 
wall and an end plate. The test cases considered for the calculation correspond to the wind tunnel 
conditions of an experiment of a NACA 0015 oscillating wing conducted in the 7 x 10 foot wind 
tunnel of the U.S. Army Aeroflightdynamics Directorate, located at the NASA Ames Research 
Center [19]. The flow conditions of the experiment are as follows: free-stream Mach Number, 
^oc = 0-29 and Reynolds number. Re = 1.95 x 10 6 based on the chord of the wing and free-stream 
speed. Four mean angles of a 0 = 4°, 11°, 15°, and 17 n are considered with an amplitude of pitch 
of a! = ±4.2° around the mean angle and a reduced frequency, k = 0.1. The instantaneous angle 
of attack of the wing is given by «(/) = ot 0 + a, sin(a>/). 

In the experiments of Piziali [19], separate two-dimensional and three-dimensional data was 
measured on a rectangular wing of aspect ratio five. Two-dimensional data was measured on the 
wing that spanned the wind tunnel wall and an end plate that replaced the tip cap. The midspan 
section at which the two-dimensional data was measured was found to be free from the cross-flow 
effects and was nearly two-dimensional. Several independent checks were made by the experimen- 
talist to assertain that this is true, including tuft survey, surface pressure measurements at several 
spanwise stations, and oil flow visualization. 

The present calculations are done for free-air conditions and, therefore, the wind tunnel walls 
are not included in the calculations. The unsteady calculations for oscillating airfoils are usually 
started from the steady state solution at the lowest angle of attack of the pitch cycle. However, 
to check the accuracy of the solution method, steady-state solutions were calculated for several 
angles of attack in a time-accurate manner. Figure 1 shows sample results of steady surface 
pressures for a = 13 J [Fig. 1(a)] and 17 [Fig. 1(b)] compared with experiments. The B- L turbulence 
model yielded satisfactory surface pressures for the case of a = 11 where the flow is essentially 



Fig. 1. Quasi-steady surface pressure distributions of NACA 0015 airfoil compared to experiments at flow 
conditions of M x ~0.29; Re = 1.95 x 10 6 ; (a) a = 13 ; (b) a = \T. 



Evaluation of turbulence models 


841 


Table l. Quasi-steady airloads 


Mean angle 

a = 13° 

a * 17° 


Airloads 

model 

C| 

Cd 

c m 

C| 

Cd 

Cm 

B-L model 

1.401 

0.0209 

0.0146 

1.624 

0.0559 

0.0267 

J-K model 

1.175 

0.0254 

0.0371 

1.093 

0.0664 

0.0391 

B-B model 

1.13 

0.0295 

0.0329 

0.7466 

0.1276 

-0.0119 

Experiments 

(Piziali) 

1.15 

0.0331 

0.028 

0.65-1.18 

0.151 

-0.028 


attached (not shown here). For the mildly separated case of a = 13° and the massively separated 
cases of a = 15 and \T\ the B-L model predicted surface pressures that had higher leading-edge 
suction peaks and milder separation. It gave higher lift and lower drag values. Therefore more 
accurate turbulence models were needed to predict satisfactory airloads. As seen in Fig. 1, the J-K 
model is able to predict the surface pressures more accurately than the B-L algebraic model at these 
conditions. Although not shown here, the B B one-equation model fared as good at 1 3 n and slightly 
better at 17 than the J-K model in predicting the overall quasi-steady airloads. The airfoil is 
completely stalled at the sc = 17° angle of attack condition. 

Table 1 lists the calculated force coefficients for these two quasi-steady flow conditions and the 
experimental measurements. The airloads are calculated by integrating the surface pressures. 
Therefore, the drag coefficient shown is that due only to pressure drag. Due to the nature of 
unsteady flow, particularly when the flow is massively separated, the calculated airloads are 
time-averaged both in experiments and computations. Therefore, the data presented tor the 
coefficients of lift (C,), drag (C d ), and pitching-moment about quarter-chord (CJ are time-aver- 
aged over a large period of time. The experimental value is also time-averaged and was read-off 
from Ref. [19]. Examination of these airloads indicates that at the lower angle of a = 13°, the J-K 
and B— B models did better in predicting the overall airloads than the B— L model, although these 
models failed to predict the drag and pitching moment accurately. At the stalled condition of 
a = 17° the flow was highly unsteady. The results are in poor agreement with each other and with 
experiment, although the agreement with experiments improved as the models got more sophisti- 
cated from B-L to B B models. 

The discussion of unsteady results is divided into three flow regimes, viz., (a) attached flow 
corresponding to a 0 = 4’; (b) light-stall case corresponding to a„= 11°; and (c) deep-stall cases 
corresponding to a 0 = 15 and 17°. In the computational procedure, the quasi-steady solution at the 
lowest angle of the pitch cycle is calculated first for each of the flow condition, and then the airfoil 
is made to execute pitching oscillations rotating about its quarter-chord point from this lowest 
angle. The evolution of unsteady flow is monitored. Typically, time periodic response of the flow 
was obtained after the second oscillation cycle. Most of the results presented in the following 
paragraphs were calculated using 361 x 71 grid with 10,000 constant time-steps per oscillating 
cycle. This number of time-steps per cycle corresponds to a nondimensional time-step of 
A/ =0.0108, based on c and a x . An explicit dissipation coefficient of c e = 0.05 was used. The grid 
size, t e , and the time-step size A / were chosen after a parametric study of these quantities using 
the J-K turbulence model. A set of values of these parameters was identified which was found to 
be optimum from the point of solution accuracy and robustness. A detailed discussion of this study 
is presented later in the section on light-stall. 

(a) Unsteady attached flow: a(f) = 4 ° + 4 . 2 ° sin(cuf) 

This flow condition essentially serves to validate the accuracy of the flow solver in calculating 
unsteady attached flow. Figure 2 presents a comparison of the calculated unsteady airloads 
obtained with different turbulence models to the experimental data. 19 In all the comparisons 
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Fig. 3. Comparison of harmonic components of unsteady pressures for different turbulence models at the 

experimental conditions of Fig. 2. 


presented in this study, the computations are compared to only unsteady loads hysteresis. A detailed 
surface pressure comparisons was not possible because of lack of availability of experimental 
surface pressure data. The experimental loads data has been averaged over 20 oscillating cycles 
but the results from calculations are not averaged but were found to be repeatable beyond the 
second cycle. The calculated loads and moments with different turbulence models compare 
favorably with each other for lift and pressure drag, but the pitching-moment appears to be 
more sensitive. The RNG, J-K, and S-A models yield relatively good comparison among 
themselves and with experimental data for lift, drag, and pitching-moment. It has been found that 
small differences in the surface pressure distributions, particularly in the trailing edge region, appear 
to produce large difference in the lift and pitching moment for the B-B model. Although, the flow 
is essentially attached for this flow condition, the B-B model shows some delay in boundary layer 
recovery in the trailing edge region during the downstroke and hence has poor comparison with 
experiments and also with loads predicted by other models. The B— L model also performs poorly 
in predicting the lift and pitching-moment for different reasons. It consistently predicts higher lift 
and lower pitching-moments compared to the rest of the models. It must be pointed out that the 
scales used in presenting the airloads in Fig. 2 are expanded to bring out the differences clearly 
for various turbulence models, but in the scales used for airloads in the rest of the study, the results 
are well within the range of experimental scatter and the differences for the various turbulence 
models. The trends of the calculated results in Fig. 2 are such that if the results are tilted slightly, 
they perhaps agree better with experiments. As observed by McCroskey and Pucci [23] in their 
experiments with NACA 0012 oscillating wing, such a trend can be attributed to wind tunnel wall 
interference effects. 

Figure 3 shows the details of unsteady pressure distributions by harmonic components, where 
the decomposition of C p is according to C p = + C p u sin(cof) + C pW cos(o>0 + C p2 sin(2cof + </>:) 

with <f> 2 as the phase-shift. In order to stretch out the leading-edge region of the airfoil, the various 
harmonics are plotted against yjx. Such a representation brings out the variations and discrepancies 
with linear theory better. Note that the differences seen in the lift and pitching-moment hysteresis 
curves for B L and B B models in Fig. 2 relatively translate into small differences in the mean and 
quadrature components in Fig. 3. No experimental pressure data was available for comparison with 
these results, but comparison of the four components w ith measurements of Ref. [23] for the NACA 
0012 oscillating airfoil for attached unsteady flow case show similar behavior for all components. 
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(b) Light-stall case: a (f)= 11 + 4.2° sin(cuf) 

Figures 4-9 show the results for the light-stall case for a 0 = 1 1 . In order to determine the 
optimum values of time-step. A/, the explicit numerical dissipation coefficient, r e , and grid size 
required for a reasonably accurate solution, parametric calculations have been performed. The J-K 
turbulence model is used for these studies. The above parameters are varied one at a time while 
keeping the others constant. For example, Fig. 4 presents results of unsteady airloads hysteresis 
for three values of At for a 361 x 71 grid and c e = 0.05 and compared to experimental data. In these 
calculations, a constant time-step At = 0.0108 is used which corresponds to using 10,000 time-steps 
per cycle. Similarly, a value of time-step = 0.0216 corresponds to 5000 time-steps per cycle. The 
maximum number of time-steps used in the present calculations is 20,000 time-steps per cycle for 
deep-stall cases. This contrasts with 50,000 to 100,000 time-steps per cycle used by other 
investigators [10, 11] in similar studies. The large number of time-steps in these investigations was 
dictated by the numerical stability of their solution method. Obviously such calculations become 
prohibitively expensive when extended to three-dimensions. As seen in Fig. 4, a time-step of 
At =0.0108 corresponding to 10,000 time-steps per cycle is a good compromise, and this is the 
number used for calculating most of the cases presented in this study. 

The numerical results presented in Fig. 5 show that the unsteady airloads and pitching-moment 
are very sensitive to the dissipation coefficient c e . The unsteady flowfields are calculated using a 
At =0.0108 on a 361 x 71 grid with J K turbulence model. For the four values of c t used, the 
upstroke results of lift, drag, and pitching-moment are nearly the same. Differences are seen for 
the downstroke. The explicit numerical dissipation apparently alters the boundary layer separation 
at the highest angle of attack around the peak of upstroke and also the reattachment process of 
boundary layer during downstroke. The evaluation and the selection of an optimal value of c c is 
done by comparing with unsteady airloads determined from experiments [19]. For example, a value 
of 0.02 for c e appears to be too small as it delays reattachment, as is evident from Fig. 5(a). Also, 
this value of c x produces nonphysical oscillatory airloads for the retreating part of the oscillating 
cycle for mean angles exceeding 11. On the other hand, a value that exceeds 0.05 seems to produce 
premature reattachment, although it appears to give better agreement for drag and pitching-mo- 
ment hysteresis during the downstroke. Choosing an optimal value of c e is therefore subjected to 
some degree of uncertainty. The uncertainties associated with the selection of e c are comparable 
to the differences produced due to variations in turbulence models. This is particularly true for the 
retreating part of the cycle. Based on these arguments, a value of c c = 0.05 is chosen as a reference 
value for the rest of the calculations. 

Figure 6 presents the unsteady airloads results calculated on different grids. Again, the J-K 
model is used along with c e = 0.05 and At =0.0108 for these calculations. The grids used in this 
study have the same wall spacing of the first grid point in the normal direction and the grid 
boundaries are located at the same distance from the airfoil. As noted before, a typical value of 
y+ for these grids is 0(1) at the midchord and there are approx. 25 points in the boundary layer. 
The stream wise grid resolution is varied for grids having 181 x 71, 361 x 71, and 671 x 71 points 
keeping the normal distribution same. The comparison of results in Fig. 6 shows sensitivity of the 
numerical solution to the grids used. A close inspection of these results shows that the 361 x 71- 
and 671 x 71 -points grids give nearly identical results. The indication of nearly grid independency 
with at least 361 points in the streamwise direction is perhaps adequate for purposes of evaluating 
turbulence models. 

Further examination of results presented in Fig. 6 indicates that the assessment of spatial grid 
resolution for the normal direction from grids 361 x 71 and 361 x 141 is less satisfactorily resolved. 
Although the grids have the same normal spacing of the first grid point at the wall as indicated 
above, its distribution in the near-wall region in the boundary layer and separated flow regions are 
quite different. This produces significant changes in c x and c m on the downstroke as seen in Fig. 
6(a) and Fig. 6(c). This suggests that although the 361-point grid in the £ direction is adequate to 
give grid-independence solution, the results from the 361 x 71 and 361 x 141 grids show that the 
71 -points grid in the t] direction is not yielding a fully grid-converged solution. In fact, the differences 
between the solutions from these two grids are no greater than the differences produced by various 
turbulence models. 
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Fig. 8. Comparison of unsteady airloads hysteresis for different turbulence models with experiments for 

the flow conditions of Fig. 4. 
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Fig. 9. Comparison of harmonic components of unsteady pressures for the J K and B B turbulence 

models at the flow conditions of Fig. 4. 


For a central-differenced numerical scheme, the explicit numerical dissipation coefficient c c has 
to be adjusted for each grid for accuracy. This is particularly true for unsteady flowfield calculations. 
In contrast, the steady flowfield simulations are not as sensitive to the range of t e used here. 
Therefore, it is not surprising that even a finer grid in the normal direction has produced poor results 
for these unsteady flows because of a single value of c e used for all the grids. It appears that from 
among the grids used here, the individual grids having 361 x 71 and 671 x 71 points are well 
matched with a dissipation coefficient of r e = 0.05 to produce acceptable results. In fact, the 
differences produced in the airloads by various grids are comparable to the difference seen for a 
given grid using different turbulence models. It may be worth mentioning that alternative to using 
a central-difference numerical scheme with added dissipation is to use a self-dissipative upwind- 
difference method [24, 25]. Such methods are found to be less sensitive to grid variation. 

A similar result of grid-sensitive study using the B-B model is presented in Fig. 7 for the same 
grids. These results are also calculated using the same values of A / and c c as those in Fig. 6. The 
unsteady airloads here show less sensitivity to the grid size. As seen, the drag and pitching-moment 
have better agreement with experiment for both upstroke and downstroke, but the lift on the 
downstroke has very poor agreement with experiment indicating that the flow reattachment is not 
complete until after the upstroke begins, which is attributable to the property of turbulence model. 
From these two grid refinement studies, the 361 x 71 grid was chosen as the optimum grid for a 
given value of At =0.0108 and c e = 0.05 and this grid is used for all further results presented. 

Using the above mentioned arguements for selecting reasonable values for A/, c e , and the 316 x 71 
points grid, unsteady flowfield solutions are calculated using the five turbulence models. Figure 8 
presents the unsteady airloads results from these solutions. Comparison of results for different 
turbulence models reveal that every model behaves differently. The chief characteristic of all these 
solutions is that all models produce trailing-edge separation. The B-L model (not shown in this 
figure) produces the least separation. The lift, drag, and pitching-moment hysteresis curves for this 
model are distinctly different from the rest of the solutions and, in particular, the lift and 
pitching-moment are in poor agreement with experiment. The lift hysteresis curves for the J-K, 
RNG, and S-A models are in good agreement with experiments, although the RNG model has 
slightly higher lift during the upstroke. In general, all models, except the B-L model, show good 
agreement with each other and with experiment for the unsteady airloads for the upstroke part 
of the cycle. As seen here, no one model has perfect agreement with experiment for all airloads. 
The B-B model shows a very slow recovery process once the boundary layer is separated and the 
reattachment is not complete until after the downstroke is complete and the upstroke begins. 
Therefore, the lift stays very low until the upstroke begins. Although not apparent in the drag curve, 
this is also seen clearly in the pitching-moment curve towards the end of downstroke. 
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The RNG and J-K models produce nearly the same extent of separation which is much larger 
than what the B-B model produces. They also appear to have similar flow recovery in the 
downstroke much better than the B-B model. Both models predict very similar drag and 
pitching-moment that are in poor agreement with experiment for the downstroke. The S-A model 
produces separation similar in extent to the B-B model, but has good flow recovery just like the 
J K model. The lift hysteresis for this model is in good agreement with experiment, like the J-K 
model, and the drag and pitching-moment are in better agreement with experiment than the J-K 
model but not as good as the B-B model. All the models except the B-B model predict poor 
pitching-moments. The B-B model produces nearly the right amount of separation and, that is the 
reason, it predicts pitching-moment correctly including its cusp-like behavior at the end of upstroke. 
The hump-like behavior at the end of downstroke is the result of poor boundary layer recovery 
due to slow reattachment as discussed before. In general, all models predict the unsteady airloads 
reasonably well for the upstroke and they all behave differently during the downstroke. The B-B 
model is the only model that predicts the pitching-moment in closer agreement with experimental 
data except for the part at the end of downstroke but nicely recovers as the upstroke begins. No 
leading edge separation initiating dynamic stall vortex was observed for this airfoil with any of the 
turbulence models. 

The NACA 0015 airfoil under investigation here, produces an unsteady flowfield separation 
initially in the trailing edge region of the upper surface of the airfoil. As the airfoil pitches up, 
during the upstroke of the oscillating cycle, the separation point moves upstream towards the 
leading edge increasing the extent of separation. The dynamic stall vortex is initiated within the 
separated flow region. Such an observation has been also made in a recent study of oscillating 
NACA 0015 airfoil by Ko and McCroskey [27]. This feature appears to be in contrast to the 
observations made for the NACA 0012 airfoil in similar oscillating environment. It is well 
documented for NACA 0012 airfoil by McCroskey et al. [23] and Chandrasekhara et al. [26] that 
the dynamic stall vortex is initiated by the leading edge separation under similar conditions. The 
phenomenon appears the same for this NACA 0012 airfoil whether the boundary layer is tripped 
in the leading edge region or not. The difference between the two airfoils is the difference in their 
leading edge curvature. Therefore, the reason for the attached flow scenario for NACA 0015 airfoil 
is the slow accelaration of flow around the leading edge region compared to NACA 0012 airfoil 
at the same flow conditions. 

There are also differences due to the compressibility for the two airfoils. It has been observed 
in the studies of Ref. [26] that shock waves are produced in the leading edge region for NACA 
0012 airfoil when the airfoil is pitched to angles exceeding 14° during the upstroke of the oscillatory 
cycle for M , greater than 0.3. These shocklets in the leading edge region trigger flow separation 
and initiate the formation of a leading edge dynamic stall vortex which will subsequently dominate 
the flow during the rest of the oscillatory cycle. In constrast, there is no evidence from either 
experiments or computations of the presence of shocks or boundary layer separation in the leading 
edge region for the NACA 0015 airfoil. The only place where separation originates for this airfoil 
is in the trailing edge region. This observation holds true even for results presented later for larger 
mean angle cases. 

Figure 9 shows the harmonic components of unsteady pressures calculated for this light-stall case 
using the J-K and B-B models. All four parts of this figure are different from those of the preceding 
attached flow case for a 0 =4° shown in Fig. 3. The large changes seen in these curves can be 
attributed to the nonlinear behavior of the flow at this mean angle of a 0 = 1 1°. The two models 
predict very similar mean and in-phase components. But the very different pitching-moments 
produced by the two models is apparent in the out-of-phase component. 

(c) Deep-stall cases: a(/) = 15° + 4.2° sin(a>t) and a(/) = 1 7° + 4.2° sin(a>/) 

In contrast to the light-stall with mild trailing-edge separation seen for ot 0 = 1 1° case, the deep-stall 
cases for a 0 = 15 and 17° are dominated by massive flow separation and highly nonlinear flow 
behavior. The boundary layer separation that originates in the trailing-edge region, during the 
upstroke of the oscillatory cycle, continues to spread upstream as the airfoil motion changes 
to downstroke. The unsteady flow behavior in this regime is characterized by the shedding of a 
large vortex-like structure during the downstroke of the cycle. This structure convects over the 
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upper surface of the airfoil and leaves the trailing edge before the downstroke part of the oscillatory 
cycle is completed. As a result, unsteady airloads far in excess of the static counterpart are generated 
during the upstroke and large amounts of hysteresis occur during the remainder of the cycle. The 
scale of the viscous-inviscid interaction zone is also large, producing a viscous layer thickness of 
the order of the airfoil chord, particularly during the vortex shedding process. 

Figure 10 shows the calculated lift, drag, and pitching-moment hysteresis loops for the a 0 = 15° 
case using different turbulence models. The results are also compared to experimental data. For 
reference, the B-L results are also shown for this case. Examination of individual curves reveals 
the following behavior. The lift hysteresis is in general agreement with calculations for all 
turbulence models except for the B-L model which shows an oscillatory-type behavior for the loads 
during downstroke. All models show very good agreement with each other and with the experiment 
for the upstroke lift curve. There are differences in the size of the dynamic stall vortex produced 
by each model and its convection downstream. Nevertheless, there are only minor differences in 
the downstroke lift curves for these models. As before, the B-B model shows clearly that it is slow 
in the recovery and reattachment process. As a result, the lift stays lower compared to the 
experimental value and the recovery is not complete until the airfoil changes its attitude from 
downstroke to upstroke motion. 

Reasonably good agreement of lift curves with experiment is not an indication to how the drag 
and pitching-moment are predicted. Unlike the light-stall results for drag and pitching-moment of 
Fig. 8, the deep-stall case shows a steep increase in drag and (negative) pitching-moment towards 
the end of upstroke of the oscillatory cycle. The J-K model which predicted good lift and drag 
curves for mean angles of a 0 = 4 and 11° also shows good agreement here for lift hysteresis with 
experiment, but predicts less than satisfactory drag and pitching-moment hysteresis curves. 
Examination of instantaneous particle flow pictures of flow for different airfoil positions during 
the downstroke reveals that the dynamic stall vortex leaves the surface much sooner for the J-K 
model than what other models predict. Except for the B-L and J- K models, all other models have 
good qualitative agreement of drag and pitching-moment hysteresis with experiment as shown in 
Fig. 10. 

The behavior of the instantaneous surface pressures during the oscillatory cycle is shown in 
Fig. 1 1 for the upper surface of the airfoil for the B-B model. As seen in this figure, the leading-edge 
suction peak continues to rise through the upstroke without stall. The peak angle of attack at the 
end of upstroke is 19.2°. This angle is about 6° above the static stall angle for this airfoil, and the 
unsteady effects extend the dynamic lift beyond the static stall angle. The leading-edge suction peak 
suddenly collapses immediately after the airfoil starts the downstroke, as revealed by the surface 
pressure distributions. Another revealing feature of this plot is that the vortex shedding phenom- 
enon manifests itself in the pressure distributions on the downstroke. 

The harmonic components of the unsteady pressures for the B-B model are presented in 
Fig. 12. Three of the four components are very different from those of the light-stall case presented 
in Fig. 9 and also of the unseparated flow presented in Fig. 3. Only the mean-component is 
qualitatively similar in shape to Fig. 9(a). The in-phase, out-of-phase, and the second harmonic 
components are all changed due to massive separation and the presence of a large-scale dynamic 
stall vortex. The growth of the second harmonic indicates the flow to be highly nonlinear. 

The various turbulence models produce different sizes of dynamic stall vortex and separated flow. 
An examination of the loci of the flow reversal point (;t s ) on the upper surface of the airfoil, 
presented in Fig. 13, shows the extent of reversed flow to vary widely at any given instant of the 
oscillatory cycle. The B-L model produces the smallest extent of reversed flow over most of the 
cycle whereas the B-B model produces the largest extent of reversed flow. The position indicating 
0° phase denotes the mean angle of oscillation. It is apparent from this figure that the B-B model 
completes the recovery process on the upstroke only when it reaches approximately the mean angle. 
For the large part of the cycle, from 15° upstroke to 15° downstroke, the RNG, B-B, and S-A 
models predict nearly the same extent of reversed flow. The massive reversed flow regions clearly 
increase the unsteady-lag effects. 

Figure 14 shows a four-part figure of the instantaneous snapshot of the flowfield during the 
oscillatory cycle computed by the B-B model. The evolution of the flow depicted clearly identifies 
the dynamic stall vortex and the extent of flow separation. The flow which separates in the 
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trailing-edge region during the upstroke cycle continues to increase in extent by moving the flow 
reversal point, S, upstream towards the leading-edge region as the airfoil continues to be pitched 
up. As seen in these snapshots, the dynamic stall vortex has peaked in its size around 17° 
downstroke and from then on it prepares to shed by moving the flow reversal point away from 
the leading edge. 

Another way of demonstrating the performance of the various turbulence models is through the 
examination of instantaneous flow pictures at any given phase in the oscillatory cycle. For example. 
Fig. 15 shows instantaneous streamline particle flow pictures for all turbulence models correspond- 
ing to the instant when the airfoil is at 16 on the downstroke. The three models, RNG, B-B, and 
S-A have very similar dynamic stall vortex structures at this instant. The B-L model has a complex 
pattern with primary and secondary vortices, whereas the J-K model has already shed the primary 
vortex at this time. All models produce multiple vortices at slightly different times during the 
downstroke. The B-L model also produces a small bubble in the leading edge region, identified 
as 51, when the airfoil has reached 14.5" in the upstroke. This bubble stays distinctly separate from 
the region containing the dynamic stall vortex and eventually merges with this at about 15° 
downstroke. The locus of the reverse flow points shown in Fig. 13 is for the point marked S in 
Fig. 15(a). 

The instantaneous streamline pictures of Fig. 15 are used only for qualitative comparison of 
different turbulence models. In fact, such a representation of flow appears to be misleading as it 
does not depict the correct picture of the unsteady flow compared to streaklines pattern. Such a 
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Fig. II. Calculated unsteady surface pressures during an oscillatory cycle for the flow conditions of 
Fig. 10 using the B B turbulence model. 
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Fig. 12. Calculated harmonic components of unsteady pressures for the flow conditions of Fig. 10 using 

the B-B turbulence model. 
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LEADING TRAILING 



Fig. 13. Locus of flow reveral point during the oscillatory cycle for different turbulence models for the 

flow conditions of Fig. 10. 
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Fig. 16. Instantaneous streakline pattern for flow at 16 downstroke using the Baldwin Barth model and 

for the flow conditions of Fig. 15. 

streakline pattern is generated by releasing particle elements at various locations in the flowfield 
and subsequently convected with local flow similar to the experimental flow visualization studies. 
Two hundred time frames of flowfield data were saved from the second oscillatory cycle and 
streaklines were constructed using the UFAT [28] program. Figure 16 shows a view of the streakline 
pattern of the flow picture shown in Fig. 15(d) at 16° downstroke for the B -B model. As seen here, 
the flow pattern and the details shown by the streaklines are phenomenally different compared to 
Fig. 15(d). The large-scale dynamic-stall vortex (VI), the pairing of vortices downstream of 
trailing-edge (V2), and a diffused pair of vortices further downstream of this (V3) is something that 
is not apparent from the instantaneous streamline patterns of Fig. 15(d). Therefore, it is necessary 
to be cautious in interpreting instantaneous steamline patterns of unsteady flowfield. 

A second case of deep-stall considered is for the mean angle of a 0 =17". This flow was calculated 
using both 10,000 and 20,000 time-steps per oscillatory cycle. It appears that for this deep-stall case 
at least 15,000 time-steps are needed to capture the important details of the flow. The results of 
unsteady airloads presented in Fig. 17 are calculated using 20,000 time-steps. Comparison of 
calculations with experiments show that all models predict the lift hysteresis fairly accurately, 
although the RNG and the J K models show oscillatory behavior during the downstroke. The 
flowfield on the airfoil for the downstroke part of the cycle is very complicated. The calculated 
results are slightly shifted from the experimental data and every model predicts separation at a 
different instant on the upstroke. Nevertheless, they all reproduce the details of the lift hysteresis 
correctly. The J-K model predicts C d and C m loops that are incorrect both for the upstroke and 
downstroke. Although the S-A model has the right trends for drag and pitching-moment, it 
produces separation too early, with the result it underpredicts the peak drag and pitching-moments. 
The B-B model calculates all the three components fairly accurately. The RNG model has good 
predictions for the upstroke but has oscillatory behavior for the downstroke. Although the 
hysteresis curves predicted by the B-B model are slightly shifted from the experimental data, it has 
demonstrated superior performance for this flow condition in predicting all unsteady airloads 
correctly. 

The B-L model is the most commonly used of all the models in computational methods and is 
the least expensive computationally. The RNG model closely follows it and is about 3% more 
expensive. The B-B model is the most expensive of the models used here and costs about 22% more 
than B-L model. The J--K and S A models are respectively 9 and 18% more expensive than B L 
model. The B-L model accounts for 12.3% of the total CPU of the numerical code. For the 
turbulence model part alone, the B-B model is about 2.5 times more expensive compared to the 
B-L model. 


CONCLUSIONS 

The unsteady, two-dimensional flowfield of an oscillating NACA 0015 airfoil is calculated using 
an implicit, finite-difference numerical method for the solution of the Navier-Stokes equations with 
an intent to evaluate the accuracy of five widely used turbulence models to calculate the unsteady 
separated flows of dynamic stall. Several unsteady flow conditions corresponding to attached flow, 
light-stall, and deep-stall cases of an oscillating wing experiment were chosen as test cases for 
calculations. First a combination of grid, time-step size, and an explicit dissipation coefficient (c c ) 
was selected based on the accuracy of solution and robustness of numerical code. The uncertainties 
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associated with the selection of c t and grid clustering in the rj direction were found to be comparable 
to the differences among various turbulence models. All calculations were performed with one set 
of these selected parameters. 

The unsteady attached flow calculated with the RNG, the J-K, and the S-A models has good 
agreement with experiments for lift, drag, and pitching-moment hysteresis. The B-L model 
performs poorly. The B-B model also predicts lift and pitching-moment hysteresis poorly for this 
attached flow case, but it has superior performance for other cases involving flow separation. 

The light- and deep-stall results are mixed. Not even one turbulence model predicts results that 
are consistent and in agreement with experiments over the three flow regimes. For the light-stall 
case, the RNG, the J-K, and the S-A models overpredict the extent of separation and therefore 
the airloads have good agreement only for the upstroke. They have a good qualitative agreement 
for lift and drag hysteresis for downstroke, but they fail to predict pitching-moments correctly. The 
B-B model predicts lower lift than experiment for the downstroke because of slow flow recovery. 
But the drag and pitching-moment are better predicted than other models. For deep-stall cases, 
the RNG, the B-B, and S -A models all predict qualitatively correct airloads hysteresis, although 
the RNG model yields oscillatory solution during the downstroke for the extreme deep-stall case. 
For the J-K model, on the other hand, the results of drag and pitching-moment hysteresis are not 
even qualitatively correct for the downstroke. 

Overall, the one-equation models provide significant improvement over the algebraic and 
half-equation models. The RNG model provides some improvement over the B-L model for nearly 
the same computational cost. Among the models considered here, the B-B model is the most 
expensive and costs about 1.22 times the B L model. The B-L model accounts for 12.3% of the 
total CPU time of the flow solver. Finally, visualization of unsteady flowfield by streaklines is more 
cumbersome but better and is consistent with the flow visualization experiments. 
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